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Abstract 

Enhanced finite elements are elements with an embedded analytical solution that can capture detailed 
local fields, enabling more efficient, mesh independent finite element analysis. The shape functions are 
determined based on the analytical model rather than prescribed. This method was applied to adhesively 
bonded joints to model joint behavior with one element through the thickness. This study demonstrates two 
methods of maintaining the fidelity of such elements during adhesive non-linearity and cracking without 
increasing the mesh needed for an accurate solution. The first method uses adaptive shape functions, where 
the shape functions are recalculated at each load step based on the softening of the adhesive. The second 
method is internal mesh adaption, where cracking of the adhesive within an element is captured by further 
discretizing the element internally to represent the partially cracked geometry. By keeping mesh adaptations 
within an element, a finer mesh can be used during the analysis without affecting the global finite element 
model mesh. Examples are shown which highlight when each method is most effective in reducing the 
number of elements needed to capture adhesive nonlinearity and cracking. These methods are validated 
against analogous finite element models utilizing cohesive zone elements. 

Introduction 

With the increasing demand for fiber reinforced composites in lightweight aerospace structures, 
adhesively bonded joints are becoming more critical than ever. Bolted and riveted joints have proven to be 
poorly suited for composite materials (Ref 1). Unlike traditional metals, brittle fibers often do not yield 
significantly to spread concentrated loads introduced by mechanical fasteners. Furthermore, bolts and rivets 
require holes in the material to be joined, which interrupts continuous fibers and introduces additional stress 
concentrations. Adhesive bonding is suitable for composite materials because it is less invasive, introduces 
load more gradually, and can often be much more cost effective. The adhesive market has indeed grown 
along with the advanced composite market, and the stmctural adhesive market in Europe has been 
forecasted to reach 67,000 tons by 2015; a growth of over 13 percent since 2008 (Ref 2). 

However, adhesively bonded joints can be problematic to model. The models often do not scale because 
of fixed thickness requirements of the adhesive layers, making individual design for each joint necessary. 
Geometric discontinuities in adherends cause stress singularities in many models, thus non-traditional 
failure criteria or evaluation methods are often required. A lack of confidence in material models, failure 
criteria, and engineering experience results in gross over-design of joints along with safety requirements 
which sometimes require secondary mechanical fasteners, jokingly referred to as “chicken bolts” by many 
engineers. 
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Adhesively bonded joints are typieally analyzed using analytieal models (elosed-form) or numerieal 
models (finite elements). Historieally, analytical models (Refs. 3 to 12) were relied upon exclusively while 
computer capability was relatively small. Analytical models are fast, and have been used to conduct 
numerous parametric studies to further the understanding of the design of joints. However, assumptions are 
often made which allow closed-form solutions but limit joint geometry and materials which can be 
accurately analyzed. Furthermore, these analytical models do not couple well with larger component and 
vehicle models, limiting their usefulness. 

Finite element (FE) models, on the other hand, are general enough to allow a wide range of geometries 
and configurations, and can even couple joint analysis with larger models. The improving speed of 
computers makes this method more viable for joint analysis. However, there are some downsides to 
modeling joints with FE models. The reentrant comers often cause a geometric singularity, and numerous 
work has been conducted to create failure theories which account for this (Refs. 13 to 19). Failure methods 
investigated include stress-based, strain based, plastic energy density, and stress measured at a characteristic 
distance from the singularity. Furthermore, the extremely thin adhesive layer limits the size of elements 
which can be used to explicitly model the adhesive. This means that there are tmly no coarse models, and 
coupling with vehicle-scale models can be problematic (Refs. 15 and 16). 

One relatively new technique for modeling progressive failure of adhesively bonded joints is 
progressive damage modeling incorporating fracture mechanics concepts. Interface elements using different 
methods such as discrete cohesive zone method (DCZM) or continuous cohesive zone method (CCZMare) 
used to resolve the stress singularity at material interfaces and reentrant geometrical comers, and allow the 
faces of the adherends to separate by treating the adhesive as a network of non-linear springs obeying a 
traction versus separation law (Refs. 20 to 26). These methods come in many different varieties, but most 
often involve stress-based initiation criterion and energy-based failure definitions. These methods have been 
shown to be extremely powerful for joints, but still have some drawbacks. Characterizing these interface 
elements requires a large amount of characterization tests, and appropriate handling of mode-mixity is also a 
subject within the cohesive zone modeling that has yet to be satisfactorily concluded (Ref 25). Additionally, 
the cohesive laws utilized require an initial, numerical fictitious stiffness to prevent separation of the plies 
before delamination initiation. Furthermore, there is a maximum element length, based on the process zone 
size, required to obtain accurate results. Cohesive zone element utilizing shape functions that are enriched 
by an analytical solution (similar to the methodology presented in this work) have been used to alleviate 
mesh dependence and size requirements (Refs. 25 and 26). Most cohesive elements are formulated 
assuming a zero-thickness interface, and thus may not be adequate to model adhesive joints, especially if the 
adhesive layer is thick. 

Another fracture based technique for modeling crack propagation that can be applied to joint analysis is 
the virtual crack closure technique (VCCT) [KRUEGER, 2004]. With VCCT fracture toughness based 
criteria are used to determine if it is energetically favorable for a crack to propagate. Propagation is 
restricted to element boundaries and typically must be known a priori. Thus, modeling joints with a finite 
thic kn ess adhesive may prove challenging with VCCT. 

All of the aforementioned methods are highly developed and have been shown to give a reasonable 
strength prediction for joints, but they are detailed models which require extremely fine meshes. Thin 
adhesive bonds, most often thinner than 1 mm, restrict the size of elements needed for the adhesive. The 
transition from the fine adhesive mesh to the coarser adherend mesh causes additional preprocessing work 
for the analyst. Therefore, joint design and analysis is typically completed after the global vehicle sizing on 
dense meshed sub-models, when design changes are expensive or impractical. 

A need exists to develop predictive tools for bonded joints that can be seamlessly coupled with large 
scale structural analyses without adding major computational demands. Such tools can be used to make 
quick mesh-independent assessments of bonded composite joints. Furthermore, they fit in into the 
computational hierarchy of virtual testing of aircraft structures (Ref 27), an area that is getting increased 
attention in the aerospace industry with the aim of lowering design cycle and certification costs. 
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A solution to this problem involves merging analytieal models with finite elements. Simplified 
struetural models ean be used to obtain shape funetions that are exaet for the assumptions of the model. 
These shape funetions ean be used to formulate stiffness matrix for the problem at hand. As long as the 
assumptions remain valid, sueh an element would give the exaet solution regardless of the number of 
elements used. 

This method has been used to ealeulate an stiffness matrix for different beam on elastie foundation 
problems (Refs. 28 and 29). More reeently, Waas and Gustafson (Ref 30) have ereated an element to 
eapture the behavior of a double overlap joint subjeeted to meehanieal and thermal loads. 

A general bonded joint finite element has been ereated (Refs. 30 to 33) wherein an entire bonded joint 
ean be modeled with a single element. This joint element eonsiders the adherends to behave like beams (or 
wide panels), and the adhesive to be made up of a bed of shear and normal nonlinear springs. The governing 
equations of this struetural model are found and solved to produee enhaneed shape funetions for the joint 
element. Furthermore, the element has been generalized to allow multiple adherend/adhesive layers and ply 
drops/thiekness tapers, providing the eapability to model a variety of joint types with very few elements. 

This model was implemented in the software Joint Element Designer, whieh was written in C# and first 
eoneeived in a joint effort between the University of Miehigan and NASA. Flowever, this method loses its 
advantage when modeling highly nonlinear adhesives and trying to eapture progressive failure. An inerease 
in elements is required for an aeeurate solution, whieh goes against the philosophy of enhaneed elements. 

This paper presents two methods whieh allow the bonded joint finite element to eapture adhesive 
non-linearities and eraeking without inereasing the mesh needed for an aeeurate solution. The first method is 
the use of adaptive shape funetions, where the shape funetions are reealeulated at eaeh load step based on 
the softening of the adhesive. The seeond method is internal mesh adaption, where eraeking of the adhesive 
within an element is represented by diseretizing a eraeked element into multiple elements in order to 
aeeurately represent the loeal, eraeked geometry. Both of these methods were implemented in the Joint 
Element Designer software. Examples are shown whieh highlight the savings in elements, eomputational 
time, and integration points needed when using these methods and when the methods are partieularly 
benefieial. The performanee of the various joint element methodologies are eompared to analogous models 
using CZM elements. 

Finally, while these methods are shown for an adhesively bonded joint but both methods have broader 
applieation. The adaptive shape funetions eould be used for any element where the material and geometrie 
properties used to obtain the shape funetions are ehanging. Updating the shape funetions within an analysis 
would improve the ability of the shape funetions to represent the deformation of the ehanging material/ 
geometry. Finally, adaptive mesh, whieh is not a new teehnique (Refs. 25, 26, 34, and 35), ean be very 
effeetive in eapturing diseontinuities within an element without ehanging the elements relationship with the 
global model. 


Formulation 

The formulation presented here is merely a summary of the formulation of the joint element, and the 
details presented are eonsidered by the authors to be required for the understanding of the eurrent study. For 
a more detailed formulation, the reader is eneouraged to eonsult prior work (Ref 36). 

Consider a strueture eonsisting of N layers of thin plates under eylindrieal bending joined together by 
N-\ thin layers of a mueh more eompliant adhesive material (Figure 1(a)). The plates are assumed to behave 
as “wide” layered eomposite Euler Bernoulli beams (henee the eylindrieal bending assumption). The axial 
displaeement of adherend /', varies linearly through the thiekness and the transverse and rotational 
displaeements, w, and ^ are assumed to be eonstant through the eross seetion. Additionally, it is assumed 
that the only signifieant stresses are the axial stresses indueed by bending and axial deformation, so other 
stress eomponents in the adherend ean be negleeted. The adhesive joining the plates is modeled as a Winkler 
foundation, where only peel and shear stress are ineluded and the stresses are assumed to be eonstant though 
the thiekness of the adhesive layer. With the displaeement behavior known in the transverse, Z;, direetion, 
the problem ean be formulated in terms of one dimension, x. 
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Figure 1 . — Adhesively bonded joint element (a) geometric parameters and (b) finite element discretization. 


Shape Functions for Linear Elastic Case 

The shape funetions used for the joint element were derived by analytieally solving the governing 
equations of the linear elastie joint model (Refs. 33 and 36). Using the prineiple of stationarity of 
potential energy and the aforementioned approximations, 2N fully eoupled governing equilibrium 
differential equations are obtained. Of the 2N governing equations, N equations eorrespond to the axial 
equilibrium, while N equations eorrespond to the transverse equilibrium. The axial displaeement 
equilibrium equations eontain seeond order derivatives, while the transverse displaeement equations have 
fourth order derivatives. The order of these equations ean be redueed and assembled into a system of first 
order constant coefficient homogeneous ordinary differential equations of the form 

u ^ = Au (1) 


where u is a vector containing the adherend centerline vectors of all of the N adherends: 


u = 







( 2 ) 


and the centerline vector of adherend i is defined as 


2 

»i=[l^iix) Wi(x) Wi(x)^^ Wi(x)^^ 


(3) 


where, x denotes the derivative with respect to x. This form of defining the centerline displacements might 
not be conventional, but it is used to lead into our solution strategy of the governing equations. Using 
state variables with higher order derivatives as is done here allows the governing equations to be reduced 
to a series of first order differential equations. 

One downside to this method worth mentioning is its dependence on solving the system of ODE’s to 
get the shape functions. This makes the extension of this method to a plate or shell-type element difficult 
because such an element would require the solution of a PDE, which can be significantly more difficult 
and time-consuming. 

The system of ordinary differential equations in Equation (1) is solved using the matrix exponential 
and rewritten in terms of the nodal degrees of freedom, q (Figure 1(b)), and shape functions, N, in the 
form: 


u=Nq 


(4) 
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The main thing to consider is that the shape functions were not prescribed, as is normally the case in finite 
elements, but derived by solving the governing equations for the stacked beam model. Therefore, the shape 
functions are a function of the geometric and material parameters of the adherend and adhesive layers: 

Material and Geometric Nonlinearities 

Since modern polymeric adhesives often display highly nonlinear material behavior, it was necessary 
to include material nonlinearities in the joint element to estimate joint strengths more correctly. A 
particularly simple nonlinear elastic stress law was chosen: 

a = a(^) ( 6 ) 

where the stress is some general function of the strain in the adhesive and adherends. Although it would 
be more correct to use an incremental flow type plasticity formulation that distinguishes loading and 
unloading stiffness, the simple nonlinear elastic relation, which assumes no permanent plastic strain, was 
chosen for several reasons. The joint element is meant to be a design tool to give general approximations, 
so it is not expected that such a tool will be used in situations requiring unloading capabilities. 
Additionally, the nature of adhesively bonded joints is such that the high stresses occur in concentrated 
form at the joint edges. Since the failing adhesive domain is eliminated in the iteration process (to be 
described later) the assumption of a nonlinear elastic type stress-strain law suffices for this modeling 
process since potential regions of “unloading” are minimal and contained in the regions which are 
eliminated. Thus, this assumption does lead to a meaningful rendition of the joint physics, yet facilitating 
an efficient (in the computational sense) solution strategy. 

The elemental stiffness matrix, k , is found by integrating the: 

k = j^B^D(q)ft/F ( 7 ) 

where B is the matrix linking the strains and nodal degrees of freedom (which contains the shape 
functions for the initial linearly elastic state: N), and D(q) is the nonlinear tangent stiffness matrix. 
Additionally, a co-rotational formulation was used to account for large rotations in the system, as 
described in References 37 to 39 and adapted to the joint element in Reference 36. 

Adaptive Shape Functions 

When the adhesive in a joint element has a nonlinear stress-strain relationship, the shape functions 
obtained for a linearly elastic adhesive may no longer be well-suited. The shape functions were derived 
based on a solution of the governing equations of the problem, so the stiffness of the materials is used to 
acquire the shape functions. As the adhesive softens, the tangent modulus changes and the shape 
functions found before loading no longer represent the solution of the governing equations based on the 
softened adhesive tangent modulus. This is illustrated in Figure 2, where a joint with a highly ductile 
adhesive is loaded and the shape functions (Figure 2(d)) are required to change as the adhesive softens. 
For highly nonlinear materials, many elements may be required to accurately capture the joint behavior. 
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Figure 2. — Demonstration of shape function progression during a loading scenario: (a) single lap joint under 
tension, (b) load versus displacement highlighting three points, and (c) tangent modulus of the adhesive and 
(d) selected shape functions at each of the three points. 


To address this issue, adaptive shape funetions were applied to the joint element. After eaeh load 
inerement, the shape funetions for the next inerement are ealeulated for a joint where the stiffness of the 
joint ehanges along the x-direetion. The funetion for the adhesive modulus, Eai(x), is based on the tangent 
stiffness of the adhesive in the prior load inerement. The peel stress in adhesive / ean be written in terms 
of the peel strain as 


( 8 ) 

and similarly for the shear stress: 

{x)r^. ( 9 ) 
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Figure 3. — To solve the governing equations for (a) an element with a softened adhesive, the joint element was 
(b) discretized into regions (in this case 6) within which the tangent modulus was considered constant. 


This causes the system of governing equations to now have non-constant coefficients (compare with 
Eq. (1)) 


u ^ = A(x)u (10) 

where the coefficient matrix, A(x) , is now a function of x. Since the coefficient matrix is non-constant, a 
semi-numerical method of solution was adopted. The domain was split into segments in which the 
coefficient matrix, A(x) , was considered constant and solved using the matrix exponential for that 
section (Figure 3). The sections were than linked together with boundary conditions, and an approximate 
solution of A(x) was found (Ref 32). 

The number of segments was determined by the number of integration points, so that the shape 
functions were refined with an increasing number of integration points. It should be acknowledged that 
the discretization involved in this method is very similar to conventional finite element discretization. 

This may serve to reduce the advantages of using the joint element, but it still allows for the joint to be 
represented with very few elements, simplifying analysis steps like mesh generation and post-processing. 

The shape functions were re-calculated at each loading step in the analysis. This adds to the 
computational time, but can be extremely advantageous in reducing the required mesh for highly 
non-linear adhesives. 


Crack Growth 

When some user-defined failure criterion is reached in some part of the adhesive layer, that portion of 
the adhesive is considered “failed” and can carry no load and has no stiffness. Setting the stress and 
stiffness of that portion of the adhesive to zero is an easy way to model the failure of the adhesive, but the 
shape functions for the joint element were not originally calculated based on a joint with failed adhesive, 
and cannot accurately model this new situation. Therefore, as with traditional finite elements, more 
elements are required to accurately find the solution where the adhesive has failed. In the case of failed 
adhesive, a great number of elements may be needed, as will be illustrated later. 

In order to increase the accuracy of the joint element after adhesive failure and crack growth, a 
method of removing the adhesive and adapting the mesh to the crack was devised. Since the joint element 
is meant to be used as a user-defined element in a larger global assembly in commercially available finite 
element software, the mesh change would have to be strictly internal to the element so that the 
surrounding model does not have to change. Therefore, a sub-assembly method was devised to handle 
adhesive failure as shown in Figure 4 and outlined in Figure 5. 
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Figure 4. — Diagram showing (a) an uncracked joint element, (b) a partially cracked element, and 
(c) a fully cracked joint element, and the addition of internal sub-elements. 



Figure 5. — Flow chart showing how cracked element sub-assembly is incorporated into joint element solution 
procedure. 


As with standard finite element analysis, a Newton-Raphson type solver is used to solve the global 
system of equations. At load increment j, the global degrees of freedom (dofs), Qj , are guessed, and this 

is broken down into elemental dofs, . After applying the prescribed displacements to the joint element, 
the adhesive is checked for failure or cracking. 

If failure in the adhesive is detected, the element is replaced by a sub-assembly with three elements as 
shown in Figure 4(b). The length of the crack, hrack, determines the lengths of the sub-assembly elements. 
The displacement of the inner nodes is unknown, therefore the sub-assembly becomes a nonlinear model 
within another nonlinear model and must be solved with its own Newton-Raphson type solution 
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procedure. A sub-assembly Newton-Raphson type solver is employed, the elemental dof are used as outer 
node boundary conditions, and the stiffness, K , and internal force vector, , for the sub-assembly 
are calculated. However, these quantities still have the inner degrees of freedom contained within. The 
force vector and stiffness matrix are then reduced using the Guyan Reduction Method (Refs. 40 to 42). 

Once the internal degrees of freedom are removed, the stiffness matrix and force vector, ky and fy”^ , can 

be considered to be that of the equivalent joint element, and are passed on to the global assembly. 

These vectors/matrices for all of the elements in the assembly are assembled, boundary conditions 
and loads are applied, and the residual, R y , (error of the initial nodal displacement guess) is calculated. If 

the residual is not within some tolerated state, a new nodal displacement “guess” is calculated based on 
the previous displacement, residual, and stiffness values and the whole cycle repeats. 

After the global system is solved, there is a check to see if the crack has grown, or if new adhesive 
failure has been detected. If this is the case, the sub-assembly is adjusted by changing the lengths of the 
sub-assembly elements, and the global system is re-solved. This is done until no new adhesive failure 
occurs and the crack is in equilibrium. A crack scaling constant, C 5 , has been introduced to speed up or 
slow down crack growth as needed, and is used in the equation 

jcur _ icur ricur _jprev \ /'I 

^ crack ~ ^ crack crack ^ crack ^ 

where is the previous crack length (prior to the global Newton-Raphson procedure) and l^rack 

the current crack length. Setting C 5 < 0 causes the crack to grow further than detected, and is useful when 
multiple iterations are needed to find crack equilibrium. Setting C 5 > 0 causes the crack to grow less than 
detected, and is necessary when crack overshoot is a concern. 

The advantage of this method is that fewer elements are needed in order to accurately capture crack 
growth. One can use the minimum elements needed to accurately capture the material and geometric 
nonlinear effects without crack growth being a factor. This can mean dramatically reducing the elements 
required, especially when there is little material nonlinearity and strains in the joint are small. 

One of the major disadvantages of this method is the increased computational time. A local nonlinear 
problem must be solved within each iteration of the global nonlinear problem. Although the local 
nonlinear problem is always limited to three elements, it can significantly increase the runtime. 
Furthermore, the global load increment is repeated if the crack grows and the sub-assemblies need to be 
created or re-meshed. Although the crack scaling parameter can significantly help in limiting the 
iterations needed to find crack equilibrium this process can still be costly. However, the costs can be 
justified if joint strength prediction is of concern. Joint strength has been identified as a controlling factor 
in the ultimate load bearing capacity of many bonded structures. 

Method 

The aim of this study was to introduce the concept of adaptive shape functions and an internal 
adaptive mesh for the joint element and show how these methods can be used to reduce the amount of 
elements or integration points needed to reach a converged solution for the failure of adhesively bonded 
joints. To demonstrate the effect that these methods have, two sets of examples will be introduced. First, 
the bonded joint element model will be compared with and without adaptive mesh and adaptive shape 
functions to show when these methods are most effective at reducing computational time and how 
effective they can be. Second, the bonded joint element model is validated against a commercial CZM 
model and the computational times will be compared. 
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Joint Element Model Comparison 


First, the bonded joint element model was used with and without adaptive shape functions (ASF) and 
adaptive mesh (AM) to see the effect of these methods. Four variations of the joint element model were 
compared (Table 1). The first one, standard joint elements, has no ASF or an AM. The shape functions 
are calculated based on the initial, undamaged state of the element. The adhesive is considered to be 
failed when the criterion is reached, and load carrying capability of the failed section of adhesive is set to 
zero in both shear and peel. The ASF model is the same as the standard model except with the shape 
functions being recalculated as the adhesive softens, and the AM model is the same as the standard model 
except that the failure of adhesive results in an introduction of internal nodes within the element and the 
position of these nodes changes as the crack in the adhesive grows. Finally, the last model, adaptive shape 
functions and mesh (ASF+M) combines both techniques. 

Two different joints are analyzed, a double cantilever beam (DCB) and a single lap joint 
(see Figure 6(a) and (b), respectively). The DCB joint was selected because the crack typically grows 
progressively through the joint and the adhesive is loaded in mode I only. The single lap joint analyzed in 
this study results in catastrophic failures (no progressive failure) and mixed mode loading. The boundary 
conditions, loading and geometry of the two joints are shown in Figure 6. 

For each joint configuration, two different adhesives are analyzed to show how effective adaptive 
shape functions and adaptive internal mesh methods are for highly plastic and brittle adhesives. The 
adhesives had an initial Young’s modulus of 4 GPa, a maximum stress of 50 MPa, an initial shear 
modulus of 1 .4 GPa, and a maximum shear stress of 29 MPa. The first adhesive displayed material 
nonlinearity only, and the stress-strain relationship was defined by a tanh function with the initial slope 
and maximum stress dictated by the aforementioned properties (see Figure 7(a)). There was no failure of 
the adhesive, it was allowed to strain indefinitely. This case is interesting because it reflects the global 
yielding failure theory introduced by Crocombe (Ref 13), which states that the upper bound of the 
strength of a joint can be found by determining the point at which the whole adhesive layer has a tangent 
modulus of zero. The second adhesive used in the study had a linear stress-strain relationship up until 
failure (Figure 7(b)). As before, the slope of the stress-strain relationship was defined previously, along 
with the stress at failure. 


TABLE 1.— FOUR JOINT ELEMENT MODEL 
VARIATIONS COMPARED 


Joint element 
model variation 

Adaptive shape 
functions 

Adaptive 

mesh 

Standard 

- 

- 

ASF 

X 

- 

AM 

- 

X 

ASF + M 

X 

X 


t A4 


I 5 mm 

t 


0.2 mm 




r|i 

1 

a) 

20 mm 

T 

50 mm 1 


j 5 mm 



1 1 i {).l mm n 


t 1 


1 

1 



1 P 

20 mm 


50 mm 



RA 


b) 

Figure 6. — Joint configurations used in the joint model comparison: (a) double 
cantilever beam (DCB) joint and (b) single lap joint (SLJ), where both joints 
are 20 mm wide (not to scale). 
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Figure 7. — ^Adhesive stress/strain relations used to demonstrate the convergence improvements using adaptive 
shape functions and adaptive mesh for adhesives displaying (a) material nonlinearity and (b) brittle material 
behavior. 


Two values were compared between the variations of the model: integration points and computational 
time to convergence. Convergence is reached when the error of the analysis is less than 1 percent. The 
error is defined in one of two ways. For joint configurations which display progressive crack growth 
(i.e., DCB joint configuration with both adhesives and single lap joint with the adhesive with nonlinear 
material properties but no failure) the error, e\, is defined as 




1 ^ |F(A,)-F^(A,)| 


( 12 ) 


where is the number of load increments, FCA, ) is the end force at load increment /', and F^(Aj) is the 
end force for the fully converged solution at load increment i. This is a measure of the average percent 
deviation of the load-displacement plot in the load-direction. When the joint fails catastrophically, or with 
no progressive crack growth, the error is taken as simply the percent deviation of the maximum load, 

or 


^2 = 


/7 _ 

max max 


100 


(13) 


where 62 is the error, is the maximum load (i.e., strength) and is the strength of the fully- 
converged solution. The analyses were all run on the same personal computer with 8 GB of RAM. While 
the actual computational time should be considered more of a qualitative property, the computational 
times of the different models relative to one another should be considered significant. 

Finally, to identify circumstances under which using adaptive shape functions becomes highly 
advantageous, a single lap joint with an adhesive containing a trapezoidal stress-strain relation and 
different levels of strain energy before failure as shown in Figure 8 was analyzed. The initial slope and 
maximum stress of the stress-strain relation is the same as that defined for the previously studied 
adhesives and the final slope was a negative value of the initial slope. The total area under the curve, 

( Wic for mode I and Wuc for mode II) was set as an integer multiple of the area the curve for the equivalent 
triangular stress-strain relation (Wai for mode I and Waii for mode II). Mode I and mode II strain energies 
were both increased, so that WJWa represents the value of WJWai and WnJWaii simultaneously. AM and 
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A 





Figure 8. — Trapezoidal stress/strain relations for (a) normal and (b) shear used to show the benefits of using 
adaptive shape functions. The length of the horizontal region of the trapezoid was extended or shortened to 
get different \NJ\Na ratios. The current energy (c) defined for mode I and II. 


ASF+M models were used, and eomputational time (i) and integration points (n) needed to reaeh 
eonvergenee was eompared in order to see how adaptive shape funetions effeets the eomputational time 
for adhesives with different amounts of eritieal strain energy. 

The adhesive was defined as failed when the following relation was satisfied: 


1 < 


W, 


IV, 


- + 


Ic 


Wn 

Wllc 


(14) 


where W/ and Wn are defined as the current strain energy, as shown in Figure 8. It should be noted that 
this is synonymous to mixed-mode laws used for cohesive zone element models, except with fracture 
energy rather than strain energy. Since the adhesive of the joint element was assigned a thickness, we 
speak of everything in terms of stress/strain rather than traction/separation. However, since the 
displacements in the adhesive are assumed to be linear in the thickness direction of the adhesive, the 
fracture energy and separation can be obtained by simply scaling the strain energy and strain by the 
adhesive thic kn ess. If the model is extended to include displacements of a higher order through the 
thickness, the traction/separation and stress/strain relations can no longer be linked. Therefore, we use the 
stress/strain relation for generality. 


Comparison With Cohesive Zone Model 

The joint element models were also compared to a cohesive zone model in Abaqus/Standard v6.1 1-1 
finite element software (Ref. 7) using pre-existing modules within the software. The implicit, static finite 
element solver was used and a finite element mesh convergence study was performed to ensure numerical 
convergence was achieved. The adherends in the Abaqus models were discretized with 2D, linear, 
quadratic, plane strain, incompatible modes, CP4EI elements. These elements were chosen to better 
represent the linear bending action in the adherends. The adherends were assigned linear elastic material 
properties. 

Cohesive elements (COH2D4) were used to represent the adhesive layers in the Abaqus models. 

These 2D elements were given a finite thic kn ess according to the thickness of the adhesive and the 
constitutive response was governed using a traction-separation description. Element tractions are related 
to the nodal separation through a traction-separation law. Typically, the area under the traction-separation 
curve is equal to the fracture toughness of the material. For these analyses, a bi-linear traction-separation 
law was chosen. A simple, maximum stress criterion was used to mark the initiation of softening within 
the element, and the power law, mixed mode implementation within Abaqus was used to couple the 
opening and shear fracture modes. Both the Abaqus and joint models used the same adhesive strengths 
used in the earlier examples. The mode I and mode II fracture toughnesses employed in the Abaqus model 
were chosen such that they were equal to the critical strain energy densities, used in the joint models, 
scaled by the thickness of the adhesive (Gjc = 625 N/m, Gnc = 601 N/mm). 
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The Abaqus built-in cohesive elements were considerably less flexible than the joint element, so the 
joint element failure criterion and constitutive relations were altered to match the cohesive elements. For 
the case of single mode loading (DCB: mode I only), the constitutive relation shown in Figure 9(a) was 
used for both models (except the cohesive zone model was scaled by the adhesive thickness to get a 
traction-separation law). The critical strain energies were the same as for the trapezoidal example when 
mymc=5. 

Since the joint element model considers the adhesive to be made of uncoupled shear and normal 
springs, the constitutive law had to be adjusted to give similar results to that of the cohesive element 
model for mixed-mode loading (SLJ). The built in cohesive zone model used a quadratic power law 
softening criterion, where softening was initiated when 



K^lcJ 


= 1 


V lie ) 


(15) 


To adapt the joint element constitutive relation to this softening criterion, the shear to peel ratio, y/ , was 
assumed to be constant throughout the loading and defined as 


y/^ 


^aMax 

I ^aMax I 


(16) 


where cr^^ax were the maximum peel and shear stresses in the adhesive found through a 

linearly elastic analysis. Finally, using Equations (15) and (16), the adjusted peel and shear stresses at 
softening, cr^ and , were found using 


f ^ 1 

1 

( ^ 
¥ 

1 

r 1 ^ 

U/C j 




\Jllc) 


and (J/c = y/Tii ^ . 


(17) 


For the cohesive zone model, the downward slopes for mode I and mode II are adjusted so that both 
modes reach zero traction when Equation (15) (peel and shear stress being replaced with mode I and 
mode II fracture energies) is satisfied, assuming that the ratio of peel to shear occurring at the onset of 
softening remains constant. For the joint element, this translates into cracking of the adhesive when: 


r 

V 


Wr 


W, 


Ic 


2 




11 


K^llcJ 


= 1 


(18) 


Assuming the ratios of peel to shear stress and peel to shear strain remain constant, this yields the 
following adjusted critical strain energies: 


[ ^ 1 


f 1 ] 

2 

+ 

f T7 2 A 
EaW 

Wlc) 


Wlc) 

[gjVuc] 


and Wjj^ 




These relation yields the mixed-mode adjusted constitutive relation shown in Figure 9(a). 


(19) 
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Figure 9. — Constitutive law joint element for comparison with cohesive elements, (a) with single-mode constitutive 

laws and (b) mixed-mode adjusted constitutive law for y/= 1.08. 

Results and Validation 

Adaptive shape funetions (ASF) and an adaptive internal mesh (AM) were introdueed to the joint 
element to reduee the number of elements needed to eapture adhesive nonlinearity and eraek growth when 
modeling an adhesively bonded joint. The following seetion illustrates the benefit of these two methods in 
redueing eomputational time, elements required, and integration points required. The first seetion 
eompares the joint element model with and without these features for different joint eonfigurations and 
adhesive types, while the final seetion eompares the model with a eommereial eode cohesive zone model. 

Joint Element Model Comparison 
Double Cantilever Beam (DCB) 

The first joint configuration to be studied was the DCB configuration. For this configuration, the 
adhesive only experiences mode I loading and crack growth occurs in a progressive manner. This analysis 
is useful in showing how the different variations of the joint element model behave under progressive 
failure and crack growth. The number of elements, number of integration points, and computational time 
for both adhesives are shown in Table 3. Corresponding load versus displacement plots are shown in 
Figure 10. Since the failure is progressive and maximum load is only one integral part of the results, the 
error measurement e\ is used and convergence occurs when the error falls below 1 percent. 

The first adhesive considered was the ductile adhesive with no failure, but softening only. Since no 
crack forms for this adhesive, adaptive mesh is not needed and only standard elements and elements with 
adaptive shape functions are compared. The adaptive shape functions were extremely effective in 
reducing the number of integration points and computational time required to get a converged solution, 
with a computational time two orders of magnitude smaller. However, even using adaptive shape 
functions, the converged solution could not be captured with merely one element, but three were required. 

The second adhesive considered was a brittle adhesive, with the stress-strain relation for mode I and 
mode II remaining linear until failure. The standard joint elements and adaptive shape function models 
both performed similarly, with similar computational times. However, the adaptive mesh method was 
pivotal in reducing the computational time significantly. Since there was no nonlinearity before failure in 
the adhesive, the adaptive shape functions did not affect the analysis except by adding to the 
computational time since the shape functions had to be recalculated at every load step. 
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TABLE 2.— COMPARISON OF NUMBER OF ELEMENTS («,,), NUMBER OF INTEGRATIONS POINTS («), 
AND COMPUTATIONAL TIME {t) FOR THE CONVERGED SOLUTION (e, < I PERCENT) FOR THE 
DCB JOINT WITH DIFFERENT MODELS AND ADHESIVES 


Model 

Ductile adhesive 

Brittle adhesive 

net 

n 

t (sec) 

»el 

n 

t (sec) 

Standard joint elements 

40 

680 

4435 

40 

680 

1988 

Adaptive shape functions (ASF) 

3 

25 

71 

1 

1025 

1498 

Adaptive mesh (AM) 




1 

65 

33 

ASF +M 

— 

.... 

— 

1 

65 

112 




Figure 10. — DCB joint configuration with load versus displacement plot for the models with (a) a ductile and 
(b) a brittle adhesive. 


Single Lap Joint (SLJ) 

The single lap joint was the next joint configuration to be considered, which is a more realistic joining 
configuration containing mode I and mode II loading in the adhesive, and catastrophic crack growth 
rather than progressive crack growth, except for the case of the nonlinear adhesive with no failure. 

When catastrophic failure occurred, it was more practical to use the error value 62 , since it deals with the 
maximum load only. The number of elements, number of integration points, and computational time for 
both adhesives are shown in Table 3. Corresponding load versus displacement plots are shown in 
Figure 1 1 . 

The first adhesive considered for the single lap joint configuration was the nonlinear adhesive with no 
failure. The maximum load is reached when the adhesive becomes fully softened in shear (the tangent 
shear modulus is nearly zero along the entire length of the adhesive). Similar to the DCB configuration, 
models with adaptive mesh were not considered because no crack growth modeling was required. As can 
be seen, using adaptive shape functions for this case made a huge difference, with two orders of 
magnitude less integration points required and three orders of magnitude faster analysis. Obviously, for 
this case, it would be highly advantageous to use adaptive shape functions. 

The second adhesive considered for the single lap joint was the brittle adhesive, which is linear until 
failure. As with the DCB configuration with the same adhesive, adaptive shape functions are not helpful 
because the adhesive has no nonlinearity, so recalculating the shape functions at every step is 
unnecessary. However, the adaptive mesh model was extremely fast for this configuration. Since the 
adhesive stress-strain relation was linear, the solution could be found in one Newton-Raphson iteration, 
with the internal mesh creation only adding a small amount of additional time. 
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TABLE 3.— COMPARISON OF NUMBER OF ELEMENTS («,,), NUMBER OF INTEGRATIONS POINTS («), 
AND COMPUTATIONAL TIME (?) FOR THE CONVERGED SOLUTION (e,/e2 < 1 PERCENT) FOR THE 
SLJ JOINT WITH DIFFERENT MODELS AND ADHESIVES 


Model 

Ductile adhesive 

Brittle adhesive 

net 

n 

t (sec) 

net 

n 

t (sec) 

Standard joint elements 

80 

1360 

35255 

40 

680 

11521 

Adaptive shape functions (ASF) 

2 

10 

71 

40 

680 

12364 

Adaptive mesh (AM) 

- 



1 

33 

2 

ASF + M 

- 



1 

33 

25 



Figure 11. — SLJ joint configuration with load versus displacement plot for the models with (a) a ductile and 
(b) a brittle adhesive. 


Trapezoidal Adhesive With Variable W 

Considering the foregoing results, one ean hypothesize that the additional eost of ealeulating the 
adaptive shape functions at every load step is beneficial for adhesives with high fracture toughness. To 
further investigate this hypothesis, the adaptive mesh model was compared with and without adaptive 
shape functions for an adhesive with a trapezoidal stress-strain relation and variable amounts critical 
strain energy W, which in this study is synonymous with the fracture toughness normalized by the 
adhesive layer thickness. The strain energy for each mode (Wic and Wuc) is normalized by the strain 
energy of an adhesive with a triangular stress-strain relation with the same critical stress and slopes. The 
computational times (t) and integration points (n) needed for a converged solution (e 2 < 1 percent) for the 
two models are compared in Figure 12. 

In Figure 12 the computational time and integration points for the converged solution of the adaptive 
mesh model is normalized by that of the adaptive shape function + adaptive mesh model. For values of 
f > 1 , the ASF+M model is faster than the AM model. The same can be said about integration points, n . 
This plot shows that for low Wc values (Wc/Wa < 3), it is not advantageous to use the adaptive shape 
functions. The two models appear to have similar results for WJWa = 3, and for Wc/Wa > 3 using adaptive 
shape functions improves both the computational time and number of integration points needed. For 
WJWa = 7, the converged solution was reached with the aid of adaptive shape functions in 1/30* the 
computational time. It is expected that the time savings would continue to increase with an increase in Wc. 
With increasing toughness of modern aerospace grade adhesives, it is expected that adaptive shape 
functions would be beneficial for many, if not most, adhesive bonding configurations. 
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Figure 12. — Comparison of the computational time (f) and integration points 
(n) needed to reach convergence (62 < 1 percent) for the joint element 
model with adaptive mesh (AM) and with both adaptive mesh and 
adaptive shape functions (ASF+M) for different stress/strain relations. 

Comparison With Cohesive Zone Model 

In order to compare the joint element model with more standard adhesive joint finite element 
techniques, a DCB and SLJ were modeled using the cohesive zone elements in Abaqus and the joint 
element models (ASF+M and AM). The coarsest finite element meshes for the Abaqus model and the 
ASF+M model are shown in Figure 13 for both joint configurations. The meshes shown represent the 
coarsest meshes that provided a converged solution (ei/e 2 < 1 percent). Both Abaqus models contained 
2800 elastic, CPE4I elements, and 100 COFI2D4 cohesive zone elements. The number of elements, 
number of integration points, and computational time for both adhesives are shown in Table 4, and the 
corresponding load versus displacement plots are shown in Figure 14. Computational time was not 
explicitly compared because the vast difference between an optimized commercial code and a flexible 
research code does not provide for a good comparison. Flowever, the number of elements and total 
number of integration points should provide a good comparison. Since both joint configurations with the 
aforementioned adhesive had significant nonlinearity and adhesive cracking, both the adaptive shape 
functions (ASF) and adaptive mesh (AM) were beneficial in providing a converged solution with a 
minimum number of elements. 

When considering the adaptive shape functions, one might conclude that with enough integration 
points, the nonlinearity of any adhesive could be captured with one element. This may be true in theory, 
but numerical stability issues caused a requirement of more elements. With very long elements and highly 
non-linear tangent stiffness (especially negative stiffness), the chances increase that the matrix 
exponential used to solve Equation (10) becomes unstable and does not converge at large values of x. 
Therefore, more elements are needed to shorten the shape functions and keep them in the stable regime. 
Future work will include applying improved matrix exponential solution procedures to correct this 
problem. 
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Figure 13. — Mesh used for the DCB (a, b) and SLJ (c, d) configurations for the 
Abaqus finite element simulations (a, c) and the ASF+M model (b, d). 


TABLE 4.— COMPARISON OF NUMBER OF ELEMENTS («,,), NUMBER OF INTEGRATIONS POINTS («), 
AND COMPUTATIONAL TIME (t) FOR THE CONVERGED SOLUTION (ei/ej < 1 PERCENT) 

FOR THE JOINTS WITH DIFFERENT MODELS 


Model 

Double cantilever beam 

Single lap joint 


n 


n 

Adaptive mesh (AM) 

320 

5440 

160 

1360 

ASF + M 

20 

340 

4 

51 

Abaqus 

2800 

11400 

2800 

11400 


a) 




Figure 14. — Comparison of load versus displacement plots of different models for (a) DCB and (b) SLJ 
configurations. 
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Summary and Conclusions 

This study deals with semi-analytical finite elements with shape functions which are determined 
based on a model rather than prescribed, specifically bonded joint finite elements. Two methods were 
demonstrated to decrease the number of elements needed for a converged answer. 

The first is adaptive shape functions (ASF). Since the shape functions are determined from a model, 
the material properties are an input parameter affecting the shape functions directly. When the material 
softens, i.e., the material properties change, the shape functions can be updated based on the new softened 
properties. This causes the shape functions to adapt to the material, making them very representative of 
the actual structural behavior and allowing a single element to remain accurate even with highly nonlinear 
materials. 

The second method for decreasing the number of elements needed was using an adaptive internal 
mesh (AM). When a portion of adhesive was cracked within a joint element, the joint element was 
re-meshed with beam elements and a joint element. This re-meshing occurred internally within the 
element, requiring no alteration of the global model. This method can capture progressive failure of an 
adhesive with very few elements, but requires a nonlinear finite element model to be solved within 
cracked element. 

Joint element models with and without the two methods were compared with each other for a double 
cantilever beam and single lap joint configuration with brittle and ductile adhesives. It was shown that the 
adaptive shape function method became more beneficial with increasing ductility of the adhesive. With 
sufficient ductility, the time spent to recalculate the shape functions at every load increment became 
insignificant when compared to the reduction in mesh and integration points. The adaptive mesh method 
was beneficial for all cases with progressive and catastrophic material cracking. 

Finally, the joint element model with both methods was compared with a standard, built-in cohesive 
zone model in Abaqus. Using both adaptive shape functions and mesh, less than 3 percent of the elements 
and integration points were needed, which is a significant gain. 

Both methods have potential for broader applications within finite elements. The adaptive shape 
functions concept can be utilized whenever the shape functions have some sort of material or geometrical 
parameter with them to improve the accuracy of shape functions during changes in materials and 
geometry. Furthermore, the adaptive internal mesh can and has been used to capture the introduction and 
propagation of discontinuities within an element without requiring a global re-meshing of the model. 
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